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O ' of the time dependent radiation transport equation is obtained using the 
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^~^' as a result of radiation interaction is fully accounted for. A tridiagonal 

,.1^ I matrix system is solved at each time step to determine the hydrodynamic 
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^ ■ the point explosion problem the self similar solutions are compared with 
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^^ , radiation energy transfer is determined. Having, thus, benchmarked the 

OO ' code, convergence of the method w.r.t. time step is studied in detail and 

^~^ I compared with the results of commonly used semi-implicit method. It is 

r~^, . shown that significant error reduction is feasible in the implicit method 

(^ ' in comparison to the semi-implicit method, though at the cost of slightly 

OO ! more CPU time. 
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1 Introduction 

Radiation transport and its interaction with matter via emission, absorption and 
scattering of radiation have a substantial effect on both the state and the mo- 
tion of materials in high temperature hydrodynamic flows occurring in inertial 
confinement fusion (ICF), strong explosions and astrophysical systems pj. For 
many applications the dynamics can be considered non-relativistic since the flow 
velocities are much less than the speed of light. In order to describe properly the 



dynamics of the radiating flow, it is necessary to solve the full time-dependent 
radiation transport equation as very short time scales (i/j ~ 1/ c or tx ~ Ap/c 
corresponding to a photon flight time over a characteristic structural length /, 
or over a photon mean free path Ap ) are to be considered [2] .Two methods 
commonly used are non-equilibrium diffusion theory [3], [5] and radiation heat 
conduction approximation [T]. The former is valid for optically thick bodies, 
where the density gradients are small and the angular distribution of photons is 
nearly isotropic. The conduction approximation is only applicable when matter 
and radiation are in local thermodynamic equilibrium, so that the radiant en- 
ergy flux is proportional to temperature gradient, and for slower hydrodynamics 
time scales [T] . Use of Eddington's factor for closing the first two moment equa- 
tions is yet another approach followed in radiation hydrodynamics . Radiative 
phenomena occur on time scales that differ by many orders of magnitude from 
those characterizing hydrodynamic flow. This leads to significant computational 
challenges in the efficient modeling of radiation hydrodynamics. 

In this paper we solve the equations of hydrodynamics and the time depen- 
dent radiation transport equation fully implicitly. The anisotropy in the angular 
distribution of photons is treated in a direct way using the discrete ordinates 
method. Finite difference analysis is used for the Lagrangian meshes to obtain 
the thermodynamic variables. The hydrodynamic evolution of the system is 
considered in a fully implicit manner by solving a tridiagonal system of equa- 
tions to obtain the velocities. The pressures and temperatures are converged 
iteratively. 

Earlier studies on the non-equilibrium radiation diffusion calculations show 
that the accuracy of the solution increases on converging the non-linearities 
within a time step and increasing benefit is obtained as the problem becomes 
more and more nonlinear and faster [B] , [7] . In this work, by iteratively converg- 
ing the thermodynamic variables, we observe a faster decrease in the L2-Error 
as compared to the commonly used semi-implicit scheme. 

The organization of the paper is as follows: In section [2] we discuss the finite 
difference scheme for solving the hydrodynamic equations followed by the solu- 
tion procedure of the radiation transport equation and their coupling together 
to result in an implicit radiation hydrodynamics code. Section [3] presents the re- 
sults obtained using this fully implicit one-dimensional radiation hydrodynamics 
code for the problems of shock propagation in aluminium and the point explo- 
sion problem. These benchmark results, thus, prove the validity of the methods. 
Next, extensive results of convergence studies w.r.t. time step are presented. 
These results using the full transport equation are new, though similar stud- 
ies have been reported earlier within approximate methods [3], [8]. Finally the 
conclusions of this paper are presented in section [H 



2 Simulation model 

2.1 Implicit finite diff'erence scheme for solving the hydro- 
dynamic equations using a Lagrangian grid. 

2.1.1 Grid structure 

For hydrodynamics calculations, the medium is divided into a number of cells 
as shown in Fig. [TJ The coordinate of the i th vertex is denoted by r.i and the 
region between the (i — 1) and i th vertices is the i th cell. The density of the i 
th grid is pi and its mass is given by 

m, ==cxp» X (r^ -rf_i) (1) 

with c — l,7r, (4/3) x n and d — 1,2,3 for planar, cylindrical and spheri- 
cal geometries respectively. Velocity of the i th vertex is denoted by Ui and 
Pi,Vi,Tion,i,Teiec,i, Eion,i and Eeieci are the total pressure, specific volume, 
temperature of ions and electrons and the specific internal energy of ions and 
electrons in the i th mesh respectively. 

2.1.2 Lagrangian step 

During a time interval At the vertexes r^ of the cells move as 

n^n + u*At (2) 

u*^{l/2){u, + u,) (3) 

where u* is the average of velocity values at the beginning and end of the 
Lagrangian step, Ui and Ui , respectively. 

2.1.3 Discretized form of the hydrodynamic equations 

In the Lagrangian formulation of hydrodynamics, the mass of each cell remains 
constant thereby enforcing mass conservation. 

The Lagrangian differential equation for the conservation of momentum is 

Here, the total pressure is the sum of the electron, ion and radiation pressures 
i.e. P = Pion + Peiec + Prad- Eq- [1] cau be discretizcd for the velocity Ui at 

1 /2 1 /2 

the end of the time step in terms of the pressures, P^ and Pili in the i th and 
i + 1 th meshes respectively, after half time step as [5] 

[Plll-P^'W 



Pi+i{n+i/2 -ri) + pi{rt -n_i/2) 



(5) 



The velocity in the i th mesh Ui is determined by the pressure in the i th 
and i + 1 th meshes and hence all the meshes are connected. Mass conservation 



equation can be used to eliminate the pressures at half time step to obtain 
an equation relating the present time step velocities in the adjacent meshes as 
follows: 

The equation describing conservation of mass is 

where p is the mass density of the medium. This equation can be rewritten in 
terms of pressure using the relation, ^ — (^)s^ — ""^Tu ^here v = J{^)s 
is the adiabatic sound speed. Therefore, Eq. [B] becomes 

- = -v^pV.u (7) 

This can be written for all the one dimensional co-ordinate systems as 

-— = -v-^p -r'^u (8) 

dt '^r" dr ^ ' 

where a = 0, 1, 2 for planar, cylindrical and spherical geometries. This equation 
can be discretized to obtain the change in total pressure along a Lagrangian 
trajectory in terms of the velocity Ui at the end of the time step [5]: 

P-J =Pi + li - PiV^ — X — - (9) 

rf_i/2 n - n-i 2 

and 

pi/2 p _L 2 1 ^^ ^ rf_^iu^+i - rfu, M 

P^il = P^+i + 9»+i - P.+i^'^+i -E. X [ — ; ; ] ^ (10) 

'i+l/2 I i+l I I ^ 

Here, qi is the quadratic Von Neumann and Richtmyer artificial viscosity in the 
i th mesh [lOj . 

_ k{p,Ax^)\ dV,^2 /11^ 

* " V, ^ dt' ^^^' 



where k (~ 3) is a dimensionless constant 

.1/ 



Using Eqs. [5] and [TU] , -P^ and -P;_(_J' in Eq. [S] are eliminated to obtain a 



tridiagonal system of equations for Ui 



where 



- AiUi+i + BiU^ - CiUi^i = Di (12) 

(13) 



Di = I -\ — — : — ^ X 



2{pAr)i ?'°Li/2(^i - ^i-i) 



(15) 



^^- 2(pAr). %r_,/,(r.-r,_i) ^^'^ 

^^ = "i ^ . . s [Pi+1 + ft+1 -Pi- qt] (17) 

(pAr), 

with 
(pAr), = Pi+i{r,+i/2 - r,) + pi{ri - r,_i/2) (18) 

The energy equations, for the ions and electrons, expressed in terms of temper- 
ature are 

p[Cv.on^^ + ^^^] = — ^^^ - ^- (19) 

and 

aTe^ec , dE,i,,dV. {Pelec + Prad)dV , 



p[C. 



Velec " 



dt dV dt' V dt 

TR{T,i,,)[ER{r, T,i,,) - B(Teiec)] + ^^e (20) 



where Eion and E'eiec are the specific internal energies and V is specific volume. 
o'B.iTeiec) IS the Rosscland opacity, Eji{r, Teiec) is the radiation energy flux and 
aii{Teiec) B(Teiec) IS radiation emission rate. Pie is the ion-electron energy 
exchange term given by 

PieiTergs/cm^/fis) = 2.704 x IQ-'^^rieiec rnon 

'-M^^Z^ -xXtlK (21) 



elec 

with ion and electron temperatures expressed in keV. Further, ^rieiec and ^nion 
are the number densities of electrons and ions, M is the mass number and Z is 
the charge of the ions. Here the Coulomb logarithm for ion-electron collision is 

m 

lnA = max{l, (23 - ln[(ne/ec)°'^ rjj^:^])} (22) 

with Teiec expressed in eV. 

The discrete form of the energy equations for ions and electrons are 



,,fc-i 

rpn,k rpn—1 _ /pn,H — l AT7-n,fc . -* je ^^'' , rn.ft — i a -r 7-n,fc\ ;^n,fc— i (')'i\ 

ion^i io7i,i \ ion^i ^i ' n k—1 'ion *^i ) I Vion^i \ "^ J 



, k ^^'>i — 1 /^r»i-_l .^.Mt- Joe ^-^1^ 



Pi 



and 



rpn,k rpn,k—l 



n,k—l^n,k—l /rpn—1 rrin,k—l\ 

ri Velec,i\ eleci eleci ) 



elec,i elec,i ' a ,T-,n,fc-l 



n,fc — 1 /pn,fe — 1 j^ pn,/c— 1 j^ cn,/c— 1\ n,fc— lAT/?^,fc 

'^m -.rj^nM _ -gn.k-l-s y^elec ^ ^rad + "elec )Pi ^^i 



^j^n,k^i>y I I I AtD"'''-^ 



+i^:' "7^r' (24) 



where 



n,k ^n.k — l 

i-,n,fc— 1 r'i ^Velec,i n,k—l^n,k—l ('ot:^ 

^i = a7 + ^Ri ^^.R^ ^^^' 

s-n.k—l /^^io7i \7i,k—l /o/?^ 



rri,fc— 1 /C-C/e(ec sn,fc— 1 

elec \ OT/ )i 

C^it' = ^ac{T:^itf (28) 






with 'n' and 'k' denoting the time step and iteration index respectively. Also, 
the constants a {— Aan/c), a^ and c denote the radiation constant, Rosseland 
opacity and the speed of light respectively. Stefan-Boltzmann law, B{Teiec) = 
acT^i^^, has been used explicitly in these equations. 

2.2 Discrete ordinates method for solving the radiation 
transport equation. 

In the Gray approximation, or one group model, the time dependent radiation 
transport equation in a stationary medium is 

Idl ^ ^ 

-—+ O.V/ + {<jr{T) + as)I{r, Q, t) = a„{T)B{T) 

c ot 

+— f i{f,n;t)dn: (29) 

where I{f,fl,t) is the radiation energy flux, due to photons moving in the di- 
rection n, at space point f and time t. Here aR{T) is the one group radiation 
opacity, which is assumed to be calculated by Rosseland weighing, at electron 
temperature T (the subscript of Teiec is dropped for convenience) . As already 
mentioned, B{T) is the radiation energy flux emitted by the medium which 
is given by the Stefan-Boltzmann law B{T) = acT^. The radiation constant 
a is ~ 137 if T is in keV and c in cm/fis. This formula for the emission 
rate follows from the local thermodynamic equilibrium (LTE) approximation, 
which is assumed in the present model. The scattering cross-section as, rep- 
resenting Thomson scattering is assumed to be isotropic and independent of 



temperature. In the Lagrangian framework the radiation transport equation for 
a planar medium is 

^P^(^) + m|^ + iMT) + 'Js)I{x, fi, t) = MT)B{T) 



/ I{x, fi/, t)d^xi. (30) 



where I{x, /i, t) is the radiation energy flux along a direction at an angle cos^^{^) 
to the X axis. The term p^(-) in this equation arises due to the Lagrange 
scheme used in solving the hydrodynamic equations. 

Backward difference formula for the time derivative gives 

^ I r ri.k — l . { \:\ — \ I 1 7-n.fc n,k—\r,n.k — \ 

+^ J ^ r''{^,l)d^,l + P—^i-\cM)-^ (31) 

Here, 'n' and 'k' denote the time step and iteration index for temperature re- 
spectively as earlier. This iteration arises because the opacity <Jr{T) and the 
radiation emission rate aii{T)B{T) are functions of the local temperature T. 
The converged spatial temperature distribution is assumed to be known for the 
hydrodynamic cycle for the previous time step. Starting with the correspond- 
ing values of (Jr{T) and B{T) , denoted by a^' and i?"'°, the radiation energy 
fluxes are obtained from the solution of the transport equation Eq. [31] . The 
method of solution, well known in neutron transport theory, is briefly discussed 
below. This is used in the electron energy equation of hydrodynamics Eq. [53] 
to obtain a new temperature distribution and corresponding values of a^ and 
B"-'^. The transport equation is again solved using these new estimates and the 
iterations are continued until the temperature distribution converges. 

Finally the transport equation can be expressed in conservation form in 
spherical geometry as 

4#(-'^"') + A[(i^Z)^] + ,r,fc = Q(,,^) (32) 

r"^ or Ofi r 



with 



<J = 'Jr'"'^ + (cAt)-^ + as (33) 






where, the second term in Eq. [32j accounts for angular redistribution of photons 
during free flight. This term arises as a result of the local coordinate system 
used to describe the direction of propagation of photons. If this term is omitted, 
Eq. [32] reduces to that for planar medium and therefore a common method of 
solution can be applied. 



In the semi- implicit method to be discussed later, the transport equation 
is solved only once per time step. Then, a slightly more accurate linearization 
[T^ can be introduced in Eqs. [21] and [32] by replacing 5"'° = i?"~^ with 
Qn,i ^ Qu ^ f^j.g^ Qj.jgj. Taylor expansion yields B" = 5"^^ + idB/dT)"{T" - 
T"~^) from which (T" — r"-i) can be eliminated using Eq. [53]. However 
this modification is not necessary in the implicit method as the iterations are 
performed for converging the temperature distribution. 

To solve Eq. [35] , it is written in the discrete angle variable as [T3] 

~2~»~v ^'") ^ (Q!m+l/2^m+l/2 ^ <^m-l/2lm-l/2) + ""-^rjj = Qm (35) 

where the indices 'n','k' on I have been supressed. Here m refers to a particular 
value of /I in the angular range [-1,1] which is divided into M directions. The 
parameter Wm is the weight attached to this direction whose value has been 
fixed according to the Gauss quadrature and q;„j±i/2 are the angular difference 
coefficients. Im and /m±i/2 are the fiuxes at the centers and the edges of the 
angular cell respectively. The angle integrated balance equation for photons is 
satisfied if the " a-coefficients" obey the condition 

'^m=l['^m+l/2lm+l/2 ^ Ckm-l/2^m-l/2] = (36) 

As photons traversing along /i = ±1 are not redistributed during the flight, the 
a-coefficients also obey the boundary conditions 

"1/2 = <^M+l/2 = (37) 

For a spatially uniform and isotropic angular flux, Eq. [35| yields the recursion 
relation 

am+l/2 = Otm-1/2 - ^m^J.m (38) 

as the flux /(r, /i) is a constant in this case. 

The finite difference version of Eq. [35| in space is derived by integrating over 
a cell of volume Vi bounded by surfaces A.1^1/2 where Vi — ^t^ Jl^iL ''^^dr = 

T'(''f+i/2 ^ ''f-i/2) ^^^ ^i±i/2 = 4'''^f±i/2' '^^^ discrete form of the transport 
equation in space and angle is thus obtained as 

/^mr. . AT 1 , 2(A,+ i/2 - A,_i/2) 

~r^l.^'i+l/2Jm,i+l/2 ~ ^i-l/2-'m.i-l/2j "1 77 

Vi UJm Vi 

'^['^m+l/2lm+l/2,i ^ Ckm-1/2 An-l/2,i] + f-^m,i = Qm,i (39) 

The cell average flux and source are given by 

Im,^ = 774^ / r^Im{r)dr (40) 

and 

1 r^i+1/2 
Qm,* = T747r/ r^Q,Jr)dr (41) 



respectively, where 'j' specifies the spatial mesh. As mentioned earlier, planar 
geometry equations are obtained if the terms involving am±i/2 are omitted and 
the replacements Vi — ri^i/2 — r^_^ii and A^j^-^ji = 1 are made. Thus, both 
geometries can be treated on the same lines using this approach. The difference 
scheme is completed by assuming that the flux varies exponentially between the 
two adjacent faces of a cell both spatially and angularly so that the centered 
flux /„i.i can be expressed as [14]: 

ir,i,i = An,j-i/2 exp [--(rj+1/2 - '"1-1/2)] (42) 

/m,j = /m.j+1/2 exp [+-(rj+i/2 - r,_i/2)] (43) 

where the radii '"i+1/2 and ri_xi2 are expressed in particle mean free paths. 
These relations show that 

^m,i — ^m,i-l/2-^m,i+l/2 (44) 

for the spatial direction. Similarly for the angular direction one gets 

^m,i — Im-l/2,ilm+l/2,i (45) 

Use of these difference schemes guarantees positivity of all the angular fluxes 
if Qyn,! are positive. The symmetry of the flux at the centre of the sphere is 
enforced by the conditions 

lM+l-mS/2 = 4i,l/2> "1=1,2, , M/2 (46) 

Dividing the spatial range into L intervals, for a vacuum boundary at ri^j^i/2, 
we have 

/™,L+i/2 = 0, m-1,2, ,M/2 (47) 

i.e, at the rightmost boundary the fluxes are zero for all directions pointing 
towards the medium. Alternately, boundary sources, if present, can also be 
specified. 

An iterative method is used to solve the transport equation to treat the 
scattering term. The radiation densities at the centre of the meshes are taken 
from the previous time step, thereby providing the source explicitly. The fluxes 
Ii/2,i for all meshes do not occur in Eq. [39] as ai/2 = 0. Then the fluxes 
/3/2,i are eliminated from this equation using the upwind scheme 1^/2,1 = h,i- 
Starting from the boundary condition, viz. Eq. [47], Eq. [39] and [44] can be 
used to determine these two fluxes for all the spatial meshes 'i'. Thereafter 
together with Eq. [55] , the fluxes for all the negative values of /i„i can be solved 
for. At the center, the reflecting boundary condition given by Eq. [46| provide 
the starting fluxes for the outward sweeps through all the spatial and angular 
meshes with positive values of ^,„ • 



This completes one space-angle sweep providing new estimates of radiation 
energy flux at the mesh centers, given by: 

^S'i ='^^mIm.i/'^UJm (48) 

in m 

where the sum extends over all directions M. The mesh-angle sweeps are re- 
peated until the scattering source distribution converges to a specified accuracy. 
The rate of radiation energy absorbed by unit mass of the material in the i th 
mesh is 

s^ = (Jji, [E^^ - B^ ]/p^ ' (49) 

which determines the coupling between radiation transport and hydrodynamics. 

2.3 Implicit radiation hydrodynamics solution method 

The sample volume is divided into 'L' meshes of equal width. The initial posi- 
tion and velocity of all the vertices are defined according to the problem under 
consideration. Also the initial pressure, temperature and internal energy of all 
the meshes are entered as input. 

For any time step, the temperature of the incident radiation is obtained by 
interpolating the data for the radiation temperature as a function of time (as in 
the case of shock propagation in aluminium sheet or an ICF pellet implosion in a 
hohlraum). All the thermodynamic parameters for this time step are initialized 
using their corresponding values in the previous time step. It is important to 
note that the velocity Ui in Eqs. [5] and [T7] and position r^ in Eq. [5] are the old 
variables and remain constant unless the pressure and temperature iterations 
for this time step converge. 

The temperature iterations begin by solving the radiation transport equation 
for all the meshes which gives the energy flowing from radiation to matter. 

The ID Lagrangian step is a leapfrog scheme where new radial velocities 
Ui arise due to acceleration by pressure gradient evaluated at half time step. 
This leads to a time implicit algorithm. The first step in the pressure iteration 
starts by solving the tridiagonal system of equations for the velocity of all the 
vertices. The sound speed is obtained from the Equation of state (EOS) of that 
material. The new velocities and positions of all the vertices are obtained which 
are used to calculate the new density and change in volume of all the meshes. 
The total pressure is obtained by adding the Von Neumann and Richtmeyer 
artificial viscosity to the ion, electron and radiation pressures and solving the 
energy equations which takes into account both the energy flow from radiation 
and the work done by (or on) the meshes due to expansion (or contraction). 
The energy equations for ions and electrons are solved using the corresponding 
material EOS which provides the pressure and the specific heat at constant vol- 
ume of the material (both ions and electrons). The hydrodynamic variables like 
the position, density, internal energy and velocity of all the meshes are updated. 
The convergence criterion for the total pressure is checked and if the relative 
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error is greater than a fixed error criterion, the iteration for pressure continues, 
i.e, it goes back to solve the tridiagonal equations to obtain the velocities, po- 
sitions, energies and so on. When the pressure converges according to the error 
criterion, the convergence for the electron temperature is checked in a similar 
manner. The maximum value of the error in electron temperature for all the 
meshes is noted and if this value exceeds the value acceptable by the error cri- 
terion, the temperature iterations continue, i.e, transport equation, tridiagonal 
system of equations for velocity, etc, are solved, until the error criterion is sat- 
isfied. Thus the method is fully implicit as the velocities of all the vertices are 
obtained by solving a set of simultaneous equations. Also both the temperature 
and pressure are converged simultaneously using the iterative method. Once 
both the pressure and temperature converge, the position of the shock front is 
obtained by noting the pressure change and the new time step is estimated as 
follows: 

The time step At is chosen so as to satisfy the Courant condition which 
demands that it is less than the time for a sound signal with velocity v to 
traverse the grid spacing Ax, ^^ < C where the reduction factor C is referred 
to as the Courant number. The stability analysis of Von Neumann introduces 
additional reduction in time step due to the material compressibility |15j . 

The above procedure is repeated up to the time we are interested in follow- 
ing the evolution of the system. The solution method described above is clearly 
depicted in the flowchart given in Fig. [5] The time step index is denoted by 
'nh' and 'dt' is the time step taken. The iteration indices for electron temper- 
ature and total pressure are expressed as 'npt' and 'npp' respectively. 'Errorl' 
and 'Error2' are the fractional errors in pressure and temperature respectively 
whereas 'etal' and 'eta2' are those acceptable by the error criterion. 

2.4 Semi-implicit method 

In the semi-implicit scheme, Eq. [5] is retained and P^ is expressed as Pj = 
{Pi + Pi)/'2. wherein Pi is the pressure at the end of the time step. Starting 
with the previous time step values for Pi , the position and velocity of each mesh 
is obtained and Pi is iteratively converged using the EOS. As the variables 
are obtained explicitly from the known values, there is no need to solve the 
tridiagonal system of equations for the velocities of all the meshes [16] . Again, 
the energy flowing to the meshes as a result of radiation interaction is obtained 
by solving the transport equation once at the start of the time step, and hence 
the iterations leading to temperature convergence are absent. 
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3 Results 

3.1 Investigation of the performance of the scheme using 
benchmark problems 

3.1.1 Shock propagation in Aluminium 

In the indirect drive inertial confinement fusion, liigli power laser beams are 
focused on the inner walls of high Z cavities or hohlraums, converting the driver 
energy to x-rays which implode the capsule. If the x-ray from the hohlraum is 
allowed to fall on an aluminium foil over a hole in the cavity, the low Z material 
absorbs the radiation and ablates generating a shock wave. Using strong shock 
wave theory, the radiation temperature in the cavity T^ can be correlated to the 
shock velocity Ug. The scaling law derived for aluminium is Tj^ ~ 0.0126u'^'^^, 
where Tr is in units of eV and Ug is in units of cm/s for a temperature range 
of 100-250 eV [If. 

For the purpose of simulation, an aluminium foil of thickness 0.6 mm and unit 
cross section is chosen. It is subdivided into 300 meshes each of width 2 x 10~^ 
cm. An initial guess value of 10~^/is is used for the time step. The equilibrium 
density of Al is 2.71 gm/cc. In the discrete ordinates method four angles are 
chosen. As the temperature attained for this test problem is somewhat low, 
the total energy equation is solved assuming that electrons and ions are at the 
same temperature (the material temperature). The Equation of State (EOS) 
and Rosscland opacity for aluminium are given by 

PV 

= eT^'V'' (50) 



kr 



7-1 



l-lj^-I^Ry-i^R (-51) 



These power law functions, of temperature and density, where e,l, fJ-, v, fJ-R and pr 
are the fitting parameters, are quite accurate in the temperature range of inter- 
est P] . 

Using the fully implicit radiation hydrodynamics code, a number of simula- 
tions are carried out for different values of time independent incident radiation 
fluxes or temperatures. Corresponding shock velocities are then determined af- 
ter the decay of initial transients. In Fig. [3l we show the comparison between 
the numerically obtained shock velocities for different radiation temperatures 
(dots) and the scaling law for aluminium (line) mentioned earlier. Good agree- 
ment is observed in the temperature range where the scaling law is valid . 

Fig. [l] shows the various thermodynamic variables like velocity, pressure, 
density and material temperature after 2.5 ns when the radiation profile shown 
in Fig. [HI is incident on the outermost mesh . This radiation temperature 
profile (Fig. [5] ) is chosen so as to achieve a nearly isentropic compression 
of the fuel pellet. The pulse is shaped in such a way that the pressure on the 
target surface gradually increases, so that the generated shock rises in strength. 
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From Fig. [D we observe that the outer meshes have ablated outwards while a 
shock wave has propagated inwards. At 2.5 ns, the shock is observed at 0.5 mm 
showing a peak in pressure and density. As the outer region has ablated, they 
move outwards with high velocities. The outermost mesh has moved to 1.2 mm. 
The meshes at the shock front move inwards showing negative velocities. Also 
the temperature profile shows that the region behind the the shock gets heated 
to about 160 eV. In Fig. [6] we plot the distance traversed by the shock front 
as a function of time for the above radiation temperature profile. The shock 
velocity changes from 3.54 to 5.46 cm/ us at 1.5 ns when the incident radiation 
temperature increases to 200 eV. 

The performance of the implicit and semi-implicit schemes are compared 
by studying the convergence properties and the CPU cost for the problem of 
shock wave propagation in aluminium. The convergence properties are examined 
by obtaining the absolute L2-Error in the respective thermodynamic variable 
profile versus the fixed time step value. The absolute L2-Error in the variable 
/ (velocity, pressure, density or temperature) is defined as 

N 

L2-Error=[5](/,-/;)2]V2 (52) 

where the data f^ constitute the exact solution for At — > 0. This exact solution 
is chosen to be the result from a run using the implicit method with a small 
time step value of lO^^^/is. The summation is taken over the values in all the 
meshes. 

Fig. [71 shows the absolute L2-Error versus the time step value for veloc- 
ity, pressure, density and temperature obtained using the implicit and semi- 
implicit radiation hydrodynamics codes respectively. The semi-implicit differ- 
encing scheme fails for time steps of 10~^/zs and higher because of the viola- 
tion of CFL criterion. For time steps ^ 10~^ fis, the errors obtained from the 
two schemes are comparable. But as the time steps are reduced, the implicit 
scheme converges very fast, i.e. the errors reduce quickly, whereas the error 
becomes nearly constant in the semi-implicit scheme because of the spatial dis- 
cretization error. An initial mesh width of 2 x lO^'' cm is chosen for the above 
convergence study, which prevents further decrease in error in the semi-implicit 
scheme. Hence a reduction in the mesh width as well as the time steps is ex- 
pected to decrease the error. In Fig. [H the i2-Error per mesh for velocity 
i-6- [Sj=i(/j ~ /j )^/number of meshes]-^/^, is plotted as a function of the time 
step by keeping the ratio of time step to mesh size i.e. At/ Ax constant at 
5 X 10~^/is/cTO. The results obtained from the implicit method using a small 
time step of At — 10~*/is and mesh width of 2 x 10^® cm is chosen as the 
exact solution. Both the implicit and semi-implicit schemes show linear conver- 
gence, though the convergence rate is faster for the implicit scheme showing its 
superiority in obtaining higher accuracies. 

Fig. [9] shows that the faster convergence in the implicit method (Fig. [71) 
is attained at the cost of slightly higher CPU time. However the cost in CPU 
seconds become comparable in the two schemes for smaller time steps. All the 
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runs in this study were done on a Pentiuni(4) computer having 1GB of RAM 
operating at 3.4 GHz. 

3.1.2 Point explosion problem 

The self similar problem of a strong point explosion was formulated and solved 
by Sedov [Hj . The problem considers a perfect gas with constant specific heats 
and density po in which a large amount of energy E is liberated at a point 
instantaneously. The shock wave propagates through the gas starting from the 
point where the energy is released. For numerical simulation, the energy E is 
assumed to be liberated in the first two meshes. The process is considered at 
a larger time t when the radius of the shock front R{t) >> ro, the radius of 
the region in which energy is released. It is also assumed that the process is 
sufficiently early so that the shock wave has not moved too far from the source. 
This ascertains that the shock strength is sufficiently large and it is possible to 
neglect the initial gas pressure Pq or counter pressure in comparison with the 
pressure behind the shock wave ^ . 

Under the above assumptions the gas motion is determined by four inde- 
pendent variables, viz, amount of energy released E, initial uniform density poi 
distance from the centre of the explosion r and the time t. The dimensionless 
quantity ^ = r/R serves as the similarity variable. The motion of the wavefront 
R{t) is governed by the relationship 

i? = ^„(l) 1/42/5 (53) 

Pq 
where ^0 is an independent variable. The propagation velocity of the shock wave 
is 

D = lf^/\^)y^R-y^ (54) 

5 po 

The parameters behind the shock front using the limiting formulas for a strong 
shock wave are 

., = :^^ (55) 

Pi = -^PoD' (56) 

7 + i 

Pi=Pa (57) 

7-1 

(7 - l)piCv 

where Cy is the specific heat at constant volume and 7 = Cp/Cy is the ratio 
of the specific heats. The distributions of velocity, pressure and density w.r.t. 
the radius are determined as functions of the dimensionless variable ^ = r/R. 
Since the motion is self-similar, the solution can be expressed in the form 

u = ui{t)uiO, P = Pi{t)P{0, P^PiPiO (59) 
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where u, P and p are new dimensionless functions. The hydrodynamic equa- 
tions, which are a system of three PDE's, are transformed into a system of 
three ordinary first-order differential equations for the three unknown functions 
M, P and p by substituting the expressions given by Eq. [SS] into the hydro- 
dynamic equations for the spherically symmetric case and transforming from 
r and t to ^. The boundary condition satisfied by the solution at the shock 
front {r = R OT^ = l)isu — P~p^l. The dimensionless parameter ^Qj 
which depends on the specific heat ratio 7 is obtained from the condition of 
conservation of energy evaluated with the solution obtained. 

Also, the distributions of velocity, pressure, density and temperature behind 
the shock front are generated numerically using the hydrodynamics code without 
taking radiation interaction into account. Ideal D-T gas of density po — 1 gm/cc 
and 7 = 1.4 is filled inside a sphere of 1 cm radius with the region divided into 
100 radial meshes each of width 0.01 cm. The initial internal energy per unit 
mass is chosen as 10^ Tergs/gm for the first two meshes and zero for all the 
other meshes. An initial time step of 10~6 ps is chosen and the thermodynamic 
variables are obtained after a time 0.2 ps. As in the case of the problem of shock 
propagation in aluminium, the total energy equation is solved assuming that 
electrons and ions are at the same temperature (the material temperature). In 
Fig. [ini we compare the distribution of the functions P/Pi,u/ui, p/pi and T/Ti 
with respect to r/R obtained exactly by solving the ODEs as explained above 
(solid lines) with the results generated from our code (dots). Good agreement 
between the numerical and theoretical results is observed. As is characteristic 
of a strong explosion, the gas density decreases extremely rapidly as we move 
away from the shock front towards the centre as seen from Fig. [101 In the 
vicinity of the front the pressure decreases as we move towards the centre by a 
factor of 2 to 3 and then remains constant whereas the velocity curve rapidly 
becomes a straight line passing through the origin. The temperatures are very 
high at the centre and decreases smoothly at the shock front. As the particles 
at the centre are heated by a strong shock, they have very high entropy and 
hence high temperatures. 

The radiation interaction effects become prominent only at higher tempera- 
tures. The Rosseland opacity of a D-T (50:50) gas in terms of the gas density 
and temperature is k^ = 0.505 p^ T"^/"^ [20], [1]. Fig. [TT] shows the profiles 
of velocity, pressure, density and temperature after 2 ps with an initial specific 
internal energy of 10^ Tergs/gm deposited in the first two meshes (i.e. total 
energy deposited is 3.351 Tergs = 3.351 x lO^^ergs ), 7 = 1.23 and 200 meshes 
each of width 0.01 cm. The shock front moves faster in the pure hydrodynamic 
case as no energy is lost to the radiation keeping the shock stronger. Also the 
effect of radiation interaction is to reduce the temperature at the centre as seen 
in Fig. [IHd). 

The convergence study for the point explosion problem shows that the L2- 
Error decreases much faster for the implicit method compared to the explicit 
scheme as observed in Fig. [12] The exact solution in this case is chosen to be 
the result from a run using the implicit method with a small time step value of 
10~^/is and an initial mesh width of 0.01 cm. As the time step is decreased. 
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the non-linearities are iteratively converged for the impHcit scheme, whereas for 
the exphcit scheme, the spatial discretization error begins to dominate thereby 
preventing any appreciable decrease in the L2-Error (as observed in the problem 
of shock propagation in aluminium foil also). In the implicit method, faster 
convergence is attained at the cost of slightly higher CPU time as shown in Fig. 

m 

4 Conclusions 

In this paper we have developed and studied the performance of fully implicit 
radiation hydrodynamics scheme as compared to the semi-implicit scheme. The 
time dependent radiation transport equation is solved and energy transfer to 
the medium is accounted exactly without invoking approximation methods. To 
validate the code, the results have been verified using the problem of shock 
propagation in Al foil in the planar geometry and the point explosion problem 
in the spherical geometry. The simulation results show good agreement with the 
theoretical solutions. The convergence properties show that the L2-Error keeps 
on decreasing on reducing the time step value for the implicit scheme. However 
for the semi-implicit scheme, the L2-Error remains fixed below a certain time 
step because of the spatial discretization error showing that the implicit scheme 
is necessary if high accuracy is to be obtained. Also the price to be paid for 
obtaining such high accuracies is only slightly higher than the semi-implicit 
scheme in terms of the CPU time. The implicit scheme gives fairly accurate 
result even for large time steps whereas the semi-implicit scheme does not work 
there because of the CFL criterion. 
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Figure 1: Grid structure. 
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Read thermodynamic parameters. Time=0, nh=0 



Time=time+dt; nh=rih+l; Read input temperature, Set parameters foi' this time step. 
■ ' ■■■ ■ npt=0 ■ -.-.....-..•.■: 



npt=rfit+l; Solve radiation transport equation. Obtain energy flowing to each mesh. 

npp=0 ■"' 



npp=npp+l ; Solve tridiagonal system of equations. Obtain new u^ Tj^ pj^ dVj. 
Solve energy equations. Obtain new Pj^ Tj^Ej^C^j for ions and electrons. 




Estimate new 
time step dt,/: 
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Figure 2: Flowchart for the Implicit ID Radiation Hydrodynamics. Here, 'nh' 
is the time step index and 'dt' is the time step taken. The iteration indices 
for electron temperature and total pressure are expressed as 'npt' and 'npp' 
respectively. 'Errorl' and 'Error2' are the fractional errors in pressure and 
temperature respectively whereas 'etal' and 'eta2' are those acceptable by the 
error criterion. 
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Figure 3: Comparison of simulation data (dots) with scaling law (line) relating 
shock velocity with the radiation temperature for aluminium. 
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Figure 4: Profiles of the tliermodynamic variables like (a) velocity, (b) pressure, 
(c) density and (d) temperature in the region behind the shock as a function of 
the distance at t = 2.5 ns. The region ahead of the shock is undisturbed and 
retain the initial values of the variables. The incident radiation temperature on 
the Al foil is shown in Figure 5. 
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Figure 5: Radiation temperature profile in the hohlraum for strong shock prop- 
agation in aluminium. 
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Figure 6: Distance traversed by the shock front vs. time graph in Al foil for the 
incident radiation temperature shown in Figure 5. The two slopes correspond 
to the two plateaus in the radiation profile. 
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Figure 7: Comparison of the L2-Error Vs. time step for the shock wave propa- 
gation problem in aluminium with mesh width Ax — 2 x 10^^ cm. Convergence 
is faster in the implicit scheme. 
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Figure 8: i2-Error/mcsh in velocity Vs. time step for the shock wave propaga- 
tion problem in aluminium with Ai/Ax = 5 x lO^^^s/cm. Convergence rate is 
higher for the implicit scheme. 
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Figure 9: CPU cost Vs. time step for the shock wave propagation problem in 
aluminium with mesh width Ax — 2 x 10~^ cm. 
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Figure 10: Comparison of the simulation data obtained in the pure hydrody- 
namic case (points) with the self similar solutions (lines) for the point explosion 
problem with specific internal energy ii^ = 10^ Tergs/gm deposited in the first 
two meshes and ^ — 1.4 . 
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Figure 11: Profiles of the thermodynamic variables with and without radiation 
interaction at 2 [is for the point explosion problem with specific internal energy 
E = 10^ Tergs/gm deposited in the first two meshes. 
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Figure 12: Comparison of the L2-Error Vs. time step for the point explosion 
problem. Implicit scheme converges faster. 
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Figure 13: CPU cost for the point explosion problem with specific internal 
energy E — 10^ Tergs/gm deposited in the first two meshes for the implicit and 
semi- implicit schemes taking radiation interaction into account. 
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